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Testing  of  RANS  turbulence  models  for  stratified 
fiows  based  on  DNS  data 


By  S.  K.  Venayagamoorthy  J.  R.  Kosefff,  J.  H.  Ferzigerf  and  L.  H.  Shihf. 


1.  Motivation  and  objectives 

Stably  stratified  flows  such  as  those  in  the  atmosphere  or  in  large  water  bodies  such  as 
the  ocean,  lakes,  estuaries  and  reservoirs  are  prevalent  in  the  natural  environment.  The 
presence  of  the  buoyancy  force  due  to  the  stratification  may  have  a substantial  effect  on 
the  flow  development  and  mixing  processes,  and  hence  influence  the  distribution  of  scalar 
substances  such  as  pollutants,  suspended  sediments,  in  the  environment. 

Today,  there  exist  numerous  turbulence  models  for  calculating  turbulent  mixing  in  the 
environment.  These  models  range  from  the  simple  eddy  viscosity  models  to  the  more 
detailed  large  eddy  simulations  (LES)  and  direct  numerical  simulations  (DNS)  (Axell 
& Liungman  2001).  However,  both  DNS  and  LES  can  be  too  computationally  expen- 
sive (and  often  too  idealized)  for  most  geophysical  and  engineering  applications.  This 
limitation  has  restricted  modelers  to  RANS  approaches  commonly  based  on  turbulent 
kinetic  energy  (TKE)  closure  schemes.  The  most  widely  used  RANS  models  today  are 
two  equation  models  which  solve  two  transport  equations  for  the  properties  of  the  turbu- 
lence from  which  the  eddy  viscosity  can  be  computed.  The  best  known  of  these  models 
is  the  k-e  model  which  requires  the  solutions  of  the  turbulent  kinetic  energy  equation  (a 
component  of  essentially  all  current  multi-equation  models)  and  dissipation  of  turbulent 
kinetic  energy  equation  (Ferziger  et  al.  2003). 

In  most  geophysical  flows,  turbulence  occurs  at  the  smallest  scales  and  one  of  the  two 
most  important  additional  physical  phenomena  to  account  for  is  stratification  (the  other 
being  rotation).  In  this  paper,  the  main  objective  is  to  investigate  proposed  changes 
to  RANS  turbulence  models  which  include  the  effects  of  stratification  more  explicitly. 
These  proposed  changes  were  developed  using  a DNS  database  on  stratified  and  sheared 
homogenous  turbulence  developed  by  Shih  et  al.  (2000)  and  are  described  more  fully  in 
Ferziger  et  al.  (2003).  The  data  generated  by  Shih  et  al.  (2000)  (hereinafter  referred  to  as 
SKFR)  are  used  to  study  the  parameters  in  the  k-e  model  as  a function  of  the  turbulent 
Froude  number,  Frk-  A modified  version  of  the  standard  k-e  model  based  on  the  local 
turbulent  Froude  number  is  proposed.  The  proposed  model  is  applied  to  a stratified  open 
channel  flow,  a test  case  that  differs  significantly  from  the  flows  from  which  the  modified 
parameters  were  derived.  The  turbulence  modeling  and  results  are  discussed  in  the  next 
two  sections  followed  by  suggestions  for  future  work. 
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2.  Turbulence  modeling  based  on  DNS  data 

2.1.  The  Data 

The  DNS  data  of  SKFR  are  used  to  develop  parameterizations  of  modeling  coefficients 
typically  found  in  RANS  models  (e.g.  C^)  as  functions  of  quantities  that  define  the  local 
state  of  the  turbulence  {Frk  etc.). 

SKFR  performed  DNS  of  stratified  homogeneous  turbulent  shear  flows.  The  data  have 
been  extensively  discussed  by  SKFR  and  hence  the  discussion  will  not  be  repeated  here. 
However,  it  is  important  to  realize  that  the  data  provide  all  of  the  properties  of  the 
turbulence  up  to  second  order  statistics.  The  subset  of  the  data  used  in  this  study  consists 
of  the  highest  initial  Taylor  microscale  Reynolds  number  runs  {Re\o  = 89).  The  value 
of  Re\o  = 89  is  relatively  high  for  direct  numerical  simulations  but  not  high  enough 
to  produce  results  that  are  independent  of  Reynolds  number  effects.  It  is  nevertheless, 
high  enough  that  the  effects  of  the  other  parameters  should  be  accurately  represented. 
All  the  data  used  for  this  study  were  taken  from  the  latter  parts  of  the  SKFR  runs  in 
order  to  ensure  that  the  turbulence  was  fully  developed.  The  physical  time  was  non- 
dimensionalized  as  St  where  S is  the  shear  rate  defined  as: 


S = 


dz 


(2.1) 


As  discussed  by  Ferziger  et  al.  (2003),  the  turbulence  does  not  become  fully  developed 
until  sometime  later  than  St  > 2.  Most  of  the  runs  were  continued  until  St  = 12-14. 
In  our  present  study,  only  the  data  for  times  between  St  = 8 and  the  end  of  the  run 
were  used.  To  render  the  plots  less  confusing,  the  data  is  averaged  over  this  time  period. 
In  total,  37  runs  were  used  to  derive  the  results  given  below.  Each  run  is  characterized 
by  the  initial  Reynolds  number  (which  is  the  same  for  all  the  runs)  and  the  gradient 
Richardson  Rig  (defined  in  equation  2.2)  which  has  a fixed  value  for  each  run. 


Rig  = 


7V2 

52 


where  N is  the  buoyancy  frequency  defined  aa: 


N = 


(2.2) 


(2.3) 


2.2.  k-e  model 

k-e  is  a commonly  used  two  equation  model,  of  which  many  variations  have  been  sug- 
gested. Here,  we  base  our  proposed  modifications  on  the  standard  version  of  the  k-e 
model.  In  the  presence  of  stratification,  the  turbulent  kinetic  energy  equation  can  be 
written  as: 

Dk 

-^  = P-e-B  + Dk  (2.4) 

where  P is  the  rate  of  production  of  TKE,  which  is  given,  for  simple  shear  flows  like  the 
SKFR  flows,  by: 

P — —u'w'S  (2.5) 

B is  the  buoyancy  flux  given  by: 


(2.6) 


RANS  Models  for  Stratified  Flows 


129 


c. 

C.i 

Ce2 

CTk 

0-£ 

0.09 

1.44 

1.92 

1.0 

1.3 

Table  1.  Values  of  constants  in  the  k-e  model  (Rodi,  1980) 


Dk  is  the  transport  term  (equal  to  zero  for  the  SKFR  flows)  modeled  using  the  gradient- 
diffusion  hypothesis  as: 


Dk 


_ d f vt  dk 
dz  \ffit  dz 


(2.7) 


and  €,  the  rate  of  dissipation  of  the  turbulent  kinetic  energy,  is  modeled  in  an  analogous 
way  to  equation  (2.4): 


Dt 


Pe 


Be 


(2.8) 


where  is  the  transport  term  (equal  to  zero  for  the  SKFR  flows)  modeled  again  using 
the  gradient-diffusion  hypothesis  as: 


The  eddy  viscosity  is  then  given  by: 


I't 


(2.9) 


(2.10) 


and  the  eddy  diffusivity  is: 


Kt  = 


I't 

Prt 


(2.11) 


where  Prt  is  the  turbulent  Prandtl  number.  The  constants  that  are  commonly  used  are 
given  in  Table  1. 

The  value  for  the  buoyancy  parameter  C^s  is  a matter  of  much  discussion.  Various 
values  have  been  suggested,  e.g.  Rodi  (1987)  suggests  a value  of  0 < Ces  < 0.29,  Baum 
& Caponi  (1992)  suggest  — 1.14  and  Burchard  & Baumert  (1995)  have  argued  that 
Ce3  should  have  negative  values.  It  turns  out  from  our  studies  that  this  is  a very  sensitive 
parameter  for  stratifled  flows. 


2.3.  Proposed  parametrizations 


Several  parameters  have  been  suggested  to  characterize  the  effects  of  the  stratification, 
the  most  obvious  of  which  is  the  gradient  Richardson  number  defined  in  equation  (2.2). 
However,  Ferziger  et  al.  (2003)  argue  that  this  is  not  the  best  choice  as  it  represents  the 
forcing  rather  than  the  properties  of  the  turbulence.  As  pointed  out  by  SKFR,  the  tur- 
bulent Proude  number  gives  better  correlations.  It  can  be  defined  based  on  the  quantities 
computed  from  the  k-e  model  as: 


(2.12) 


Further,  it  is  also  noted  that  there  is  a gradient  Richardson  number  at  which  turbulence 
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energy  neither  grows  nor  decays  called  the  stationary  Richardson  number  Rigs-  Holt  et 
al.  (1992)  showed  that  Riga  is  a function  of  the  Reynolds  number  at  least  at  smaller 
values  of  Re.  If,  in  equations  (2.4)  and  (2.8),  we  set  Dk/Dt  = De/Dt  = 0,  and  eliminate 
e,  then  for  homogeneous  flows,  the  following  relation  is  obtained: 


(2.13) 


where  Rif  a is  the  flux  Richardson  number.  Applying  the  gradient  diffusion  modeling 
concept  to  both  B and  P in  equation  (2.13)  yields: 


Rif  a — 


Riga 

Tn 


(2.14) 


where  it  is  found  that  Prt  w 1 in  the  stationary  flow  cases.  Hence  we  assume  that  Ri fa 
and  Riga  are  equivalent  in  the  discussion  that  follows. 

We  choose  to  insist  that  be  zero  at  the  stationary  state  in  order  to  satisfy  the 
constraint  that  the  model  predicts  the  existence  of  a stationary  state  under  the  conditions 
determined  by  SKFR.  Rearranging  equation  (2.13)  and  using  the  generally  accepted  value 
of  0.25  for  the  stationary  Richardson  number  at  high  turbulence  Re3molds  number,  we 
get: 


1 Rif  a 


(2.15) 


The  functional  dependence  of  the  stationary  Richardson  number  Riga  on  the  turbulence 
Reynolds  number  Re^  defined  in  equation  (2.16)  was  given  by  SKFR  as  shown  in  equation 
(2.17). 


€1/ 


(2.16) 


Rif,  « Riga  - ^ ^ (2.17) 

Any  of  the  parameters  in  the  model  e equation  can  be  allowed  to  vary  as  a function 
of  the  turbulent  Proude  number  but  only  one  parameter  (or  combination  of  parameters) 
can  be  derived  from  the  model  dissipation  equation  and  the  DNS  data,  enforcing  a 
limitation  on  the  choices  that  can  be  made.  In  this  study,  the  parameters  that  are  chosen 
to  depend  on  the  stratification  are  <7^3  and  while  Cd  and  0^2  are  independent  of 
the  stratification.  We  do  however  include  the  effect  of  the  Reynolds  number  on  as 
discussed  earlier  through  equations  (2.15)  and  (2.17). 

The  buoyancy  parameter  can  be  calculated  from  DNS  data  using  the  model  e 
equation  as  a function  of  Fr*,.  To  accomplish  this,  we  compute  [l/e)de/dt  from  the  SKFR 
data  by  fltting  a least  square  straight  line  to  log  e as  a function  of  time  over  the  time 
range  used  in  all  of  the  data  fitting.  We  then  solve  equation  (2.8)  for  using  equation 
(2.15)  for  Ca2  (as  discussed  by  Ferziger  et  al,  2003).  The  resulting  values  are  plotted  in 
Figure  1.  Unfortunately,  there  is  no  clear  trend  in  the  data  to  suggest  any  definitive  fit. 
However,  the  data  indicates  that  is  of  order  unity  in  the  strongly  stratified  region 
and  possibly  shows  a slight  increase  with  increasing  Frk-  The  fit  we  used  for  this  study 
is  also  shown  in  Figure  1 and  is  given  by  equation  (2.18).  This  correlation  is  based  on 
the  observation  that  the  mixing  efflcency  peaks  at  about  Frk  of  0.4  - 0.5  for  the  SKFR 
data  and  should  be  zero  there.  It  should  be  noted  that  the  data  does  not  extend  into 
the  very  strongly  stratified  regime  (i.e.  very  low  Frk  values).  It  is  a regime  that  has  a 
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Figure  1.  Buoyancy  parameter  Ce3  as  a function  of  the  turbulent  Proude  number,  o:  values 
calculated  from  the  SKFR  data, ; fit  given  by  (2.18). 


mixture  of  internal  waves  and  turbulence  and  care  has  to  be  exercised  in  applying  any 
turbulence  model  in  this  regime. 


1.44 

1.44  - 9.6(Prfc  - 0.35) 
6A{Frk  - 0.5) 

1.92 


for  Frk  < 0.35 
for  0.35  < Frk  < 0.50 
for  0.5  < Frk  < 0-80 
for  Frk  > 0.80 


(2.18) 


We  note  that  Figure  1 indicates  that  a better  fit  to  the  data  could  be  given  but  we  felt 
it  important  to  repeat  the  condition  that  be  zero  at  Frk  = 0.5. 

The  eddy  viscosity  parameter  obtained  from  the  data  is  plotted  as  function  of 
Frk  in  Figure  2.  The  value  of  obtained  from  the  SKFR  data  at  large  Frk  (weak 
stratification)  is  lower  than  the  typical  value  of  0.09.  However,  we  fit  the  data  such  that 
reaches  the  asymptotic  value  of  0.09  at  high  Frk  values.  Thus: 

{0.125Fr^  + 0.014Frfc  for  Fr*,  < 0.35 

0.006(Fr*  - 0.35)7(0.02  + 0.1(Frfe  - 0.35))  + 0.02  for  0.35  < Frk  < 0.60 
0.08tanh(Frfc)  + 0.01  for  Frk  > 0-60 

(2.19) 

The  scalar  transport  can  be  modeled  using  the  turbulent  Prandtl  number  defined  in 
equation  (2.11).  We  plotted  the  turbulent  Prandtl  number  as  a function  of  the  turbulent 
Proude  number  as  shown  in  Figure  3.  The  curve  fit  shown  in  the  figure  is: 


f 1.4  for  Frk  < 0.35 

\ 1.4  — 0.55(1  — exp(— 7(Frfc  — 0.35))  for  Frk  > 0.35 


(2.20) 


where  we  have  chosen  to  keep  the  Prandtl  number  constant  in  the  highly  stratified 
regions. 
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Figure  2.  The  eddy  viscosity  coefficient  as  a function  of  the  turbulent  Proude  number,  o: 
values  derived  from  the  SKFR  data, : correlation  (2.19) 


Figure  3.  The  turbulent  Prandtl  number  as  a function  of  the  turbulent  Proude  number,  o: 
values  derived  from  the  SKFR  data, : correlation  (2.20) 


3.  Results  and  discussion 

In  order  to  test  the  proposed  parameterization  and  compare  it  with  other  models,  it 
is  essential  that  we  apply  it  to  a flow  that  is  significantly  different  from  the  homogenous 
flows  from  which  the  parameters  were  derived.  We  chose  the  stratified  open  channel  flow 
as  our  test  case.  This  is  a flow  for  which  direct  numerical  simulations  and  large  eddy 
simulations  have  been  performed  by  Garg  et  al.  (2000)  and  Shih  (2003). 

The  tests  of  the  model  were  performed  using  a 1-D  water  column  model  called  GOTM 
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(General  Ocean  Turbulence  Model)  developed  by  Burchard  et  al.  (1999).  The  various 
turbulence  closure  schemes  in  GOTM  are  based  on  a modular  format  that  enables  easy 
incorporation  in  3-D  ocean  circulation  models  and  also  allows  for  refinements  or  exten- 
sions to  the  turbulence  models.  The  source  code  is  available  to  the  public  on  the  Internet 
website  www.gotm.net.  We  incorporated  our  proposed  model  into  the  GOTM  code.  Since 
none  of  the  built  in  stability  functions  in  GOTM  modify  their  parameters  based  on  depth, 
we  had  to  adapt  the  code  to  modify  the  parameters  in  k-e  model  discussed  above  to  be 
depth  dependent  based  on  the  local  turbulent  Proude  number. 

The  test  case  we  use  is  somewhat  artificial.  It  is  a pressure-gradient  driven  open  channel 
flow  in  which  the  density  is  held  fixed  at  both  the  lower  solid  boundary  and  the  upper 
free  surface  and  in  some  ways  similar  to  the  experiments  done  by  Komori  et  al.  (1983). 
The  Reynolds  number  Rer  and  the  Richardson  number  Rir  based  on  the  shear  velocity 
for  all  the  test  runs  were  682  and  31  respectively.  Previous  studies  by  Garg  et  al.  (2000) 
have  shown  that  LES  produces  results  that  are  in  good  agreement  with  DNS  results  for 
open  channel  flows.  Hence  we  use  an  LES  run  with  identical  conditions  to  those  in  our 
test  case  to  assess  the  predictions  of  our  proposed  model  and  that  of  the  standard  k-e 
model  with  constant  stability  functions. 

The  flow  was  first  allowed  to  develop  to  a converged  solution  (without  stratification 
eff'ects)  using  the  standard  k-e  model  after  which  the  stratification  was  imposed.  The 
velocity  and  density  profiles  at  the  initial  state  (i.e.  after  spin  up)  are  shown  in  Figure 
4.  Also  shown  in  Figure  4,  are  the  turbulent  kinetic  energy  and  dissipation  profiles  to- 
gether with  the  turbulent  Froude  number  and  the  turbulent  Reynolds  number  profiles. 
Superimposed  on  these  figures  are  the  instantaneous  LES  profiles. 

A total  of  four  runs  were  done  using  different  combinations  of  the  modifications  (out- 
lined above)  to  the  k-e  model  so  as  to  determine  the  most  suitable  combination  of  param- 
eters that  can  match  as  closely  as  possible  the  LES  results.  The  test  runs  discussed  in  this 
paper  are  outlined  in  Table  2.  Run  1 is  based  the  standard  k-e  model  with  = 1.44, 
while  run  2 uses  the  proposed  parametrizations.  In  run  3,  we  vary  only  the  eddy  viscosity 
parameter  and  turbulent  Prandtl  number  Prt  and  keep  both  Ce2  and  constant 
as  in  run  1.  Run  4 is  an  arbitrary  case  that  was  chosen  to  investigate  the  effects  of  just 
varjdng  the  turbulent  Prandtl  number  using  the  correlation  given  by  equation  (3.1). 

The  velocity,  density,  turbulent  kinetic  energy,  dissipation  of  the  turbulent  kinetic  en- 
ergy, turbulent  Proude  number  and  turbulent  Reynolds  number  profiles  at  non-dimensional 
times  tu^/h  = 4.1  and  tu^/h  = 14.6  for  the  standard  k-e  model  (run  1)  and  the  full 
modified  version  (run  2)  together  with  the  LES  profiles  are  shown  in  Figures  5 and  6 
respectively. 

The  full  modified  version  of  the  proposed  model  (denoted  by  run  2)  underperforms 
the  k-e  model  in  terms  of  predicting  the  velocity  profile  in  the  channel  especially  at  later 
times,  as  seen  in  Figme  6,  even  though  it  appears  to  capture  the  free  surface  velocity 
better.  However,  it  does  significantly  better  in  predicting  the  evolution  of  the  density 
profile,  especially  during  the  transient  stages  of  the  flow. 

In  run  3,  both  C(2  and  Cj3  are  held  at  constant  values  equal  to  the  ones  used  for  the 
k-e  test  run  1.  It  can  be  seen  from  the  density  profile  in  Figure  7 that  the  mixing  is 
now  predicted  better  for  the  highly  stratified  regions  (i.e.  low  Proude  numbers)  in  this 
run.  However,  this  case  underperforms  the  standard  k-e  case  for  most  of  the  lower  half 
of  the  channel  depth.  The  discrepancy  between  the  LES  and  the  RANS  results  given  by 
runs  1 and  3 for  the  lower  half  of  the  channel  depth  indicates  that  in  general  the  k-e 
model  allows  too  much  mixing  to  occur  from  the  bottom  boundary.  There  is  a strong 
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Figure  4.  Initial  profiles  (a)  velocity;  (b)  density;  (c)  tke;  (d)  dissipation;  (e)  turbulent  Proude 
number  Fr*;  (f)  turbulent  Reynolds  number  Rek-  • — : LES  results, : k-e  results. 


Run 

C.2 

Ce3 

Pn 

Remarks 

1 

0.09 

1.92 

1.44 

0.85 

Standard  k-e  model 

2 

eq.(2.19) 

eq.  (2.15) 

eq.  (2.18) 

eq.  (2.20) 

Pull  modified  version 

3 

eq.(2.19) 

1.92 

1.44 

eq.  (2.20) 

Varying  and  Prt  only 

4 

0.09 

1.92 

1.44 

eq.  (3.1) 

Varying  Prt  only 

Table  2.  Test  runs  used  for  evaluating  the  predictions  by  the  modified  k-e  model 
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Figure  5.  Profiles  at  tu,/h  — 4.1:  (a)  velocity;  (b)  density;  (c)  tke;  (d)  dissipation;  (e)  Frk\ 
(f)  Rck-  — : LES  results; : run  1;  : run  2. 

possibility  that  the  prescribed  asymptotic  value  of  0.85  that  the  SKFR  data  suggests 
for  the  turbulent  Prandtl  number  is  low  even  though  there  are  good  engineering  data 
that  suggest  this  value  is  reasonable.  However,  other  observations  such  as  those  due  to 
Hogstrom  (1996)  suggest  a value  of  Prt  = 1.  In  run  4,  we  use  a different  correlation  for 
Prt  such  that  it  asymptotes  to  1.0  at  high  Frk  (see  equation  3.1).  The  results  as  shown 
in  Figure  7 indicate  closer  agreement  with  the  LES  results  compared  to  runs  1 and  3 
respectively.  . 

Prt  = 0.4exp(— 2.5Frfe)  + 1.0  for  all  Frk  (3.1) 


4.  Conclusions  and  futinre  work 

Modifications  of  the  k-e  model  to  account  for  stratification  based  on  the  SKFR  data 
have  been  suggested  and  tested  in  this  study.  The  test  runs  done  on  the  stratified  open 


Figure  6.  Profiles  at  tu,/h  = 14.6:  (a)  velocity;  (b)  density;  (c)  tke;  (d)  dissipation;  (e) 

(f)  Rck.  — : LES  results; : run  1;  •••  ; run  2. 

channel  flow  highlight  the  importance  of  correctly  modeling  the  turbulent  Prandtl  number 
Prt  as  well  as  the  buoyancy  parameter  C^s.  The  results  suggest  that  the  turbulent  Prandtl 
number  should  be  close  to  unity  for  the  neutrally  stable  flows.  Further,  it  appears  that 
the  buoyancy  parameter  has  to  be  prescribed  as  a value  of  the  order  of  C^i  in  order 
to  correctly  model  the  effects  of  the  buoyancy  force  in  the  dissipation  equation.  The 
lack  of  a clear  trend  in  as  a function  of  Fr*  from  the  DNS  data  suggests  that  the 
strategy  of  setting  0^3  equal  to  a constant  is  probably  the  most  reasonable  approach.  The 
results  of  this  paper  are  useful  only  in  the  regimes  of  weak  to  moderate  stratification. 
Clearly  further  test  cases  should  be  performed  where  the  effects  of  the  stratification  are 
more  pronounced  and  where  Frk  is  small  over  larger  fraction  of  the  flow  depth.  This  will 
allow  us  to  evaluate  the  effects  of  varying  (7^  as  function  of  Fr/.  more  precisely.  Further 
work  is  also  required  to  capture  the  boundary  layer  dynamics  properly  so  that  momentum 
balance  can  be  achieved  properly.  We  are  also  unable  to  comment  on  the  transport  models 


Figure  7.  Velocity  profiles  (left)  and  density  profiles  (right).  — : LES  results; : runl;  •••  ; 

run  2;  - • - : run  3;  - • - in  bold  : run  4 
at  tu,/h  = 4.1  (top)  and  tu^/h  — 14.6  (bottom). 


(which  may  also  depend  on  stratification)  as  the  SKFR  database  is  for  homogeneous  flows. 
It  is  important  to  obtain  data  for  strongly  stratified  and  inhomogeneous  flows  and  apply 
them  to  the  development  of  better  models  for  stably  stratified  flows. 
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